#ifndef _MCMCHISTOGRAM_H
#define _MCMCHISTOGRAM_H

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <assert.h>
#include "Vector.h."

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
// histogram
// This class is rewriten for the conservation of memory space.

class Histogram
{
public:
	// construction & destructions
	Histogram();
	~Histogram();
	void Free(void);

	// operators
	Histogram&	operator=(const Histogram &a);
	Histogram&	operator+(Histogram &a);
	double		operator*(Histogram &a);
	FILE & operator>>(FILE &strm);
	FILE & operator<<(FILE &strm);

	virtual	void WriteToTxtFile(FILE &strm);
	virtual	void ReadFromTxtFile(FILE &strm);

	// operations
	virtual void	Set(const double x1, const double x2, const double d1);
	inline	void	SetCutOff(const bool bcut){cut_off=bcut;}
	inline	void	SetBin(const double x, const double p);
	inline	double	GetBin(const double x);
	inline	void	AddBin(const double x, const double p);
	inline	void	IncrementBin(const double a){AddBin(a,1.0);};
	inline	void	SetIndexedBin(const int i, const double p);
	inline	void	AddIndexedBin(const int i, const double p);
	virtual void	Normalize(void);
	inline	double	GetValueIndex(const int i){return x0+d*i;};

	virtual void	Init(const double a);
	double			Entropy(void);
	double			DistTV(Histogram &a);
	double			DistKL(Histogram &a);
	double			DistToGaussian(void);
	double			ComputeMean(void);
	inline	double	GetMeanBinNum(void){return mean;};
	inline	double	GetMeanDomain(void){return (double)(mean*d+x0);};
	mcmcVector<double>	BinDifference(Histogram &a);
	inline  int		GetIndex(const double x);
	virtual int	TotalCodingBits(void);

private:
	bool			AllocateSpace(const int size);

public:
	double			*pbins;					// data
	int			count;					// count of pdf
	double			x0;						// starting value
	double			d;						// interval
	double			outof_lower_bound;		// the number of values that are out of the lower bound
	double			outof_upper_bound;		// the number of values that are out of the high bound

protected:
	double			mean;
	bool			cut_off;				// true: values that are out of bounds are considered
};


inline int Histogram::GetIndex(const double x)
{
	return (int)floor((x-x0)/d);
}

inline void Histogram::SetBin(const double x, const double p)
{
	// get index
	int i = GetIndex(x);

	if (i>=0 && i<count)
		pbins[i] = p;
	else if (!cut_off)
	{
		if (i < 0)
			outof_lower_bound = p;
		else
			outof_upper_bound = p;
	}
}

inline double Histogram::GetBin(const double x)
{
	// get index
	int i = GetIndex(x);

	if (i>=0 && i<count)
		return pbins[i];
	else if (!cut_off)
	{
		if (i < 0)
			return outof_lower_bound;
		else
			return outof_upper_bound;
	}
	else
		return 0.0;
}

inline void Histogram::AddBin(const double x, const double p)
{
	// get index
	int i = GetIndex(x);
	if (i>=0 && i<count)
		pbins[i] = pbins[i]+p;
	else if (!cut_off)
	{
		if (i < 0)
			outof_lower_bound += p;
		else
			outof_upper_bound += p;
	}
}

inline	void Histogram::SetIndexedBin(const int i, const double p)
{
	assert(i>=0 && i<count);
	pbins[i] = p;
}

inline	void Histogram::AddIndexedBin(const int i, const double p)
{
	assert(i>=0 && i<count);
	pbins[i] += p;
}

//////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////
// Gibbs parameters
class Gibbs_B: public Histogram
{
public:
	// constructions & destructions
	Gibbs_B(){;};
	~Gibbs_B(){;};

	// operators
	Gibbs_B& operator=(Gibbs_B& a);
	friend FILE & operator>>(Gibbs_B &g, FILE &strm);
	friend FILE & operator<<(Gibbs_B &g, FILE &strm);

	// operations
	inline	void SetB4OutLowerBound(const double b){outof_lower_bound=b;};
	inline  void SetB4OutUpperBound(const double b){outof_upper_bound=b;};
	double	Projection(Histogram &hist);
	double	ComputeZ(void);
	inline	double GetBin(const double x);
};

inline double Gibbs_B::GetBin(const double x)
{
	int i=GetIndex(x);
	if (i<0)
		return outof_lower_bound;
	else if (i>=count)
		return outof_upper_bound;
	else
		return pbins[i];
}


#endif
